% This script will carry out error analysis
R = 0.5;    % 0.5Ohm
L = 0.0015; % 1.5mH

data_points = 10000;

% Go on for a time constant
time_constant = L/R;
step = time_constant/data_points;

T=150e-6; % 150us
f = 1/T;  
w_c = 2*pi*f;
A = 6;

Vin = @(t) A*cos(w_c*t);
current_initial=0;

% e^m is the integratingn factor
m = 0.5/0.0015;

% Solution is made by multiplying by integrating factor and
% then integrating both sides
current_exact = @(t) 3/(m^2+w_c^2)*(2*m*cos(w_c*t) + 2*w_c*sin(w_c*t));
Vout_exact = @(t) Vin(t) - current_exact(t)*R;

[time_array, Vout_array] = ralston(R, L, Vin, current_initial, step, data_points*step);

for j=1:data_points
	error_array(j) = Vout_exact(time_array(j)) - Vout_array(j);
	% Vout_plot(j) = Vout_exact(time_array(j));
end

% plot(time_array, Vout_array);
% figure;
% plot(time_array, Vout_plot);
% 
% max(Vout_plot-Vout_array)

loglog(time_array, error_array);