Here’s a trick to quickly convert a continuous, strictly proper transfer function into pseudocode for a discrete time implementation.  Let’s start with a first order low pass filter:

$$\frac{Y(s)}{U(s)}=\frac{1}{\tau s+1}$$

The first step is to cross multiply the terms (for simplicity I’ll omit the Y(s) and U(s) for the rest of this post).

$$\tau sY+Y=U\to \tau \dot{Y}+Y=U$$

Then solve the transfer function for the highest ordered derivative of $Y$ (meaning the $Y$ term with the most $\cdot$ above it), then integrate both sides of the equation until $Y$ (with no dots) becomes the left hand side.  The integration is with respect to time, but I’ve omitted that for clarity in the expression.

$$\dot{Y}=\frac{U-Y}{\tau }\to \int \dot{Y}=\int \frac{U-Y}{\tau }\to Y=\int \frac{U-Y}{\tau }$$

Now replace the $=\int \left(\right)$ with a $+=dt\left(\right)$ to achieve the following pseudocode

$$Y+=dt\left(\frac{U-Y}{\tau }\right)$$

Here, $dt$ is the time between filter updates, and $Y$ is a static variable whose value persists between filter updates.  What we did is solve for the derivative of the filter output ($\dot{Y}$) and then integrated the derivative using Euler’s method to get the filter’s software implementation.

We can do the same thing with a second order filter:

$$\frac{Y}{U}=\frac{{\omega }_n^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2}$$

$$\ddot{Y}+2\zeta {\omega }_n\dot{Y}+{\omega }_n^{2}Y={\omega }_n^{2}U\to \ddot{Y}=-2\zeta {\omega }_n\dot{Y}+{\omega }_n^2(U-Y)$$

Here we’ve again solved for the highest ordered derivative of $Y$.  Now what we’ll do is integrate both sides of the equation until we are again left with only $Y$ on the left hand side.

$$\int \ddot{Y}=\int -2\zeta {\omega }_n\dot{Y}+{\omega }_n^2(U-Y)\to \dot{Y}=-2\zeta {\omega }_{n}Y+\int {\omega }_n^2(U-Y)$$

$$\int \dot{Y}=\int (-2\zeta {\omega }_{n}Y+\int {\omega }_n^2(U-Y))\to Y=\int (-2\zeta {\omega }_{n}Y+\int {\omega }_n^2(U-Y))$$

Let’s break the previous expression into two parts so only one $\int ()$ appears on the right hand of any equals sign.

$$k=\int {\omega }_n^2(U-Y)$$

$$Y=\int (-2\zeta {\omega }_{n}Y+k)$$

Now we pull the same trick we did before by replacing each $=\int ()$ with $+=dt()$ to arrive at the software pseudocode.

$$k+=dt({\omega }_n^2(U-Y))$$

$$Y+=dt(-2\zeta {\omega }_{n}Y+k)$$

Now we have two static variables, $k$ and $Y$, that each represent a separate state variable of the second order filter.

As an exercise for the reader (oh man, did I really just say that?), show that the second order band pass filter

$$\frac{Y}{U}=\frac{2\zeta {\omega }_{n}s}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2}$$

has the following pseudocode implementation

$$k+=dt(-{\omega }_n^{2}Y)$$

$$Y+=dt(k+2\zeta {\omega }_n(U–Y))$$