work in progress solution for exercise 1 of coursework

1 files changed, 25 insertions, 0 deletions

diff --git a/coursework17/heun.m b/coursework17/heun.m
new file mode 100644
index 0000000..ab7edfd
--- /dev/null
+++ b/coursework17/heun.m
@@ -0,0 +1,25 @@
+function [time_array, Vout_array] = heun(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/2*(grad + feval(dcurrent_dt, time_array(j)+step, current_array(j)+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
+