function [time_array, Vout_array] = ralston(R, L, Vin, current_initial, step, t_final)
	% Determine how many times to step forward
	N = round(t_final/step);

	% Initialise some arrays for the output
	current_array=zeros(1,N);
	time_array=zeros(1,N);
	current_array(1) = current_initial;
	Vout_array(1) = Vin(0) - R*current_initial;

	% Given the input voltage as a function of time we can express 
	% dI/dt (dcurrent_dt) as a function of time as well (given R and L)
	dcurrent_dt = @(t, I) (Vin(t) - I*R)/L;

	% Loop as many times as step allows and:
	% calculate the gradient at point
	% Use it to emtimate next two points

	for j=1:N-1 % loop N-1 times
		grad = feval(dcurrent_dt, time_array(j), current_array(j));
		current_array(j+1)=current_array(j) + step/4*(grad + 3*feval(dcurrent_dt, time_array(j)+2/3*step, current_array(j)+2/3*step*grad)) ; 
		time_array(j+1) = time_array(j) + step;
		Vout_array(j+1) = Vin(time_array(j+1)) - R*current_array(j+1);
	end