Attempting to use the trapezoidal rule to form a difference equation representing a circuit I have a differential equation that has been proven to be correct. 

The transfer function obtained by Laplace domain analysis and Matlab freqs match up and all is well.

The problem is somewhere in my attempt to digitize the differential equation. I am required to do this using the trapezium rule to approximate the analog signal. 

The goal is to get a digitized equation in the form of..

Vo^n+1 = alpha1(Vi^n+1) + alpha2(Vi^n) + alpha3(Vo^n)

Where each alpha value is made up of a different combination of T, R1, R2 and
C all numerical parameters.

I am following some guidelines set by my lecturer but if the approach seems wrong please let me know.

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Current incomplete solution..





Then I am trying to convert that equation to the Z domain and therefore use the freqz tool in Matlab. Please see the working out below.



Here is the Matlab plot that shows what has gone wrong. The reason for the multiple lines are varying values of capacitance for C.



The blue lines represent the validated analogue freqs(laplace domain transfer function) frequency response.

The green lines show my attempt at using the freqz(z domain discrete transfer function) frequency response.

If it is of any further interest, the code to generate those plots can be found below as well.

Any help would be sincerely appreciated, clutching at straws and the occasional hair at the moment!

Cheers

MATLAB code:

%% Comparing the impulse amlitude response using two different methods
%
%  Basically checking if freqs matches freqz

clc; clear all;

% Circuit Parameters
R1 = 2.2e3;
C = 0.01e-6;
R2 = 10e3;

fs = 10000;   % Sampling frequency (Hz)
T = 1/fs;     % Timestep

% Coefficients
den = (T*R1 + T*R2 + 2*R1*R2*C);

N = 1000;

% Digital parameters
y = zeros(N,1);
y(1) = 1;
result = zeros(N,1);

% Analogue parameters
fa = 0:fs/N:fs-fs/N;
w = 2*pi*fa;
j =sqrt(-1);


lastValue = 0;

clf;
figure(1);

hold on;

for i=1:5

% Analogue Recalculations
b = 1/(R1*C);
a = (R1 + R2)/(R1*R2*C);
Ha = freqs(b,[1 a], w);

% Digital Recalculations
alpha1 = (T*R2)/den;
alpha2 = alpha1;
alpha3 = (2*R1*R2*C - T*R1 - T*R2)/den;

[Hd,fd] = freqz([alpha1 alpha1],[1 alpha3],N,fs);

for n=1:N-1

 currentValue = alpha1*y(n+1) + alpha2*y(n) - alpha3*lastValue;
 result(n) = currentValue;
 lastValue = currentValue; 
end

C = 1.5*C;

plot(fa,abs(Ha),fd,abs(Hd)),xlabel('Frequency'),ylabel('H(s)');

end

hold off;
The Laplace transform of the original differential equation is
$$Cs V_o(s) = V_i(s) \left( \frac{1}{R_1} \right) - V_o(s) \left( \frac{1}{R_1} + \frac{1}{R_2} \right)$$
$$R_1 R_2 C s V_o(s) = R_2 V_i(s) - (R_1 + R_2) V_o(s).$$
In the Laplace domain, the ratio $V_o(s)/V_i(s)$ is the transfer function for the circuit,
$$H(s) = \frac{ V_o(s) }{ V_i(s) } = \frac{ R_2 }{ R_1 + R_2 + R_1 R_1 C s }
= \frac{K}{1 + R_p C s},$$
where $K = R_2 / (R_1 + R_2)$ is the voltage divider ratio for $R_1$ and $R_2$
in a series connection, and $R_p = R_1 R_2 / (R_1 + R_2)$ is the equivalent
resistant for $R_1$ and $R_2$ in a parallel connection.
We can obtain the $z$-domain transfer function for an equivalent digital filter
using the bilinear transformation,
$$s = \frac{2}{T} \left( \frac{z-1}{z+1} \right)$$.
In this transformation, the constant is set to $2/T$, which will produce a digital filter that is an approximation to the analog filter that is equivalent to approximating the convolution integral by the trapezoidal rule.  I assume this is what you are looking for.
By making this substitution, the transfer function for the digital filter is
$$H(z) = \frac{K}{1 + R_p C \left( \dfrac{2}{T} \right) 
\left( \dfrac{z-1}{z+1} \right)}
= \dfrac{ KT (z+1) }{ T(z+1) + 2 R_p C (z-1) }$$.
This can be further simplified to
$$H(z) = \frac{V_o(z)}{V_i(z)} =
\frac{ \left( \dfrac{KT}{T + 2 R_p C} \right) \left( 1 - z^{-1} \right) }
{1 + \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) z^{-1}}.$$
Taking the inverse $z$-transform of this equation, we obtain a
difference equation for the digital filter,
$$Y(z) \left[ 1 + \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) z^{-1} \right]
= X(z) \left[ \left( \dfrac{KT}{T + 2 R_p C} \right) \left( 1 - z^{-1} \right) \right].$$
$$y[n] = \left( \dfrac{KT}{T + 2 R_p C} \right) \left( x[n] - x[n-1] \right)
- \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) y[n-1].$$
You should be able to use this equation to recursively compute $y[n]$.