1 %{
 2 RK4
 3 Octave 5.1.0
 4 John Bryan 2019
 5 %}
 6
 7
 8 clear all;
 9 close all;
10
11
12 function PlotFunction(w_array,h,tn,i)
13    figure (i);
14    k=0:h:tn;
15    exact=zeros(length(k));
16    exact=19*exp(0.85*k)
17    plot(k, exact, '.g+');
18    hold on
19    plot(k, w_array,'.ob');
20    handle = legend ("exact", ["RK4 h=" num2str(h)]);
21    legend (handle, "location", "northwest");
22    legend boxoff;
23    set(gcf,'InvertHardCopy','off');
24    filename=["rk4_plot" int2str(i) ".png"];
25    print (gcf, "-S500,500", filename, '-dpng'); 
26    return
27 endfunction
28
29
30 function v = fo(t,y) 
31    v=.85*y;
32    return; 
33 endfunction
34
35
36 function  w_array = rungekutta(f,w,tn,h) 
37    t=0;
38    N=tn/h
39    w_array=zeros(1,N+1);
40    w_array(1)=w;
41    for i=[1:N]
42       k1=h*f(t,w);
43       k2=h*f(t+h/2,w+k1/2);
44       k3=h*f(t+h/2,w+k2/2);
45       k4=h*f(t+h,w+k3);
46       w=w+(k1+2*k2+2*k3+k4)/6;
47       t=t+h;
48       w_array(i+1)=w;
49    endfor
50    return ;
51 endfunction
52
53
54 f=@fo;
55 y0=19;
56 w=y0;
57 tn=3;
58 h=[0.5,0.25,0.1];
59 for i=[1:1:length(h)]
60    w_array=rungekutta(f,w,tn,h(i))
61    PlotFunction(w_array,h(i),tn,i)
62 endfor