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