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

data_points = 100000;
h_divisions = 1000;

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

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

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


A = 6*R/(R^2 + (w*L)^2);
B = 6*w*L/(R^2 + (w*L)^2);

current_exact = @(t) A*cos(w*t) + B*sin(w*t) - A*exp((-R/L)*t);
Vout_exact = @(t) Vin(t) - current_exact(t)*R;

for k=1:h_divisions

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

	for j=1:data_points/k
		error_array(j) = Vout_exact(time_array(j)) - Vout_array(j);
    end
	max_array(k) = max(error_array);
	max_t(k) = k*step;
end

loglog(max_t, max_array);
title('Order of error')
xlabel('Time step size')
ylabel('Maximum error')