program ABVSS_Example_2; {------------------------------------------------------------------------------} { } { The program solves the initial value problem } { y' = 0.5y, y(0) = 1, } { with the exact solution y = exp(0.5t), by interval methods of } { Adams-Bashforth type with finding a step size for each integration step. } { Note: In order to speed up calculations one can use iadd, isub, imul and } { idiv functions from our IntervalArithmetic32and64 unit (see } { ABVSS_Example_2_ASMD program). In this program we use overloading } { operators also implemented in our unit. } { } {------------------------------------------------------------------------------} {$APPTYPE CONSOLE} {$R *.res} uses System.SysUtils, System.DateUtils, IntervalArithmetic32and64; function modify_time (const time1, time2 : TDateTime) : string; var code, days, hour_numbers : Integer; hours, time, tstr : string; h, m, s, ms : Word; begin DecodeTime (time1-time2, h, m, s, ms); Str (h, tstr); if h<10 then time:='0'+tstr else time:=tstr; Str (m, tstr); if m<10 then time:=time+':0'+tstr else time:=time+':'+tstr; Str (s, tstr); if s<10 then time:=time+':0'+tstr else time:=time+':'+tstr; Str (ms, tstr); if ms<10 then time:=time+'.00'+tstr else if ms<100 then time:=time+'.0'+tstr else time:=time+'.'+tstr; days:=DaysBetween(time1,time2); if days>0 then begin hours:=Copy(time,1,2); Val (hours, hour_numbers, code); hour_numbers:=24*days+hour_numbers; Str (hour_numbers, hours); Result:=hours+Copy(time,3,10) end else Result:=time end {modify_time}; function F (const T, Y : interval) : interval; begin Result:=Y/2 end {F}; function PSI1(const T, Y : interval) : interval; begin Result:=Y/4 end {PSI1}; function PSI2(const T, Y : interval) : interval; begin Result:=Y/8 end {PSI2}; function PSI3(const T, Y : interval) : interval; begin Result:=Y/16 end {PSI3}; function PSI4(const T, Y : interval) : interval; begin Result:=Y/32 end {PSI4}; procedure IntervalAB1 (const T1, Y1, ih, delta_t, delta_y : interval; out Y : interval); var PSI, Th : interval; begin Th.a:=0; Th.b:=ih.b; PSI:=Y1+Th*F(delta_t,delta_y); PSI:=PSI1(T1+Th,PSI); Y:=Y1+ih*F(T1,Y1); Y:=Y+(ih*ih/2)*PSI end {IntervalAB1}; procedure IntervalAB2 (const T1, Y1, T2, Y2, ih1, ih, delta_t, delta_y : interval; out Y : interval); var PSI, Th, g : interval; begin Th.a:=-ih1.b; Th.b:=ih.b; PSI:=Y1+Th*F(delta_t,delta_y); PSI:=PSI2(T1+Th,PSI); Y:=F(T1,Y1)-F(T2,Y2); Y:=F(T1,Y1)+(ih/ih1)*Y/2; Y:=Y1+ih*Y; g:=(1/3+(ih1/ih)/2)/2; Y:=Y+ih*ih*ih*g*PSI end {IntervalAB2}; procedure IntervalAB3 (const T1, Y1, T2, Y2, T3, Y3, ih2, ih1, ih, delta_t, delta_y : interval; out Y : interval); var PSI, Th, g : interval; begin Th:=ih1+ih2; Th.a:=-Th.b; Th.b:=ih.b; PSI:=Y1+Th*F(delta_t,delta_y); PSI:=PSI3(T1+Th,PSI); Y:=F(T2,Y2)-F(T3,Y3); Y:=(ih1/ih2)*Y; Y:=F(T1,Y1)-F(T2,Y2)-Y; g:=(ih/ih1)*((ih+ih1)/(ih1+ih2)); g:=(1-ih/(3*(ih+ih1)))*g; Y:=g*Y/2; g:=F(T1,Y1)-F(T2,Y2); g:=((ih/ih1)*g)/2; Y:=F(T1,Y1)+g+Y; Y:=Y1+ih*Y; g:=(ih1/ih)*((ih1+ih2)/ih)/2; g:=((2*ih1+ih2)/ih)/3+g; g:=(1/4+g)/6; Y:=Y+ih*ih*ih*ih*g*PSI end {IntervalAB3}; procedure IntervalAB4 (const T1, Y1, T2, Y2, T3, Y3, T4, Y4, ih3, ih2, ih1, ih, delta_t, delta_y : interval; out Y : interval); var PSI, Th, g : interval; begin Th:=ih3+ih2+ih1; Th.a:=-Th.b; Th.b:=ih.b; PSI:=Y1+Th*F(delta_t,delta_y); PSI:=PSI4(T1+Th,PSI); Y:=F(T3,Y3)-F(T4,Y4); Y:=(ih2/ih3)*Y; Y:=F(T2,Y2)-F(T3,Y3)-Y; g:=(ih1/ih2)*((ih1+ih2)/(ih2+ih3))*Y; Y:=(ih1/ih2)*(F(T2,Y2)-F(T3,Y3)); Y:=F(T1,Y1)-F(T2,Y2)-Y-g; g:=(ih+ih1+ih2)/(ih1+ih2+ih3); g:=((ih+ih1)/(ih1+ih2))*g; Y:=(ih/ih1)*g*Y; g:=ih/(ih+ih1+ih2); g:=((1-ih/(2*(ih+ih1)))/6)*g; g:=(1-(ih/(3*(ih+ih1))))/2-g; Y:=g*Y; g:=(ih1/ih2)*(F(T2,Y2)-F(T3,Y3)); g:=F(T1,Y1)-F(T2,Y2)-g; g:=(ih/ih1)*((ih+ih1)/(ih1+ih2))*g; g:=(1-ih/(3*(ih+ih1)))*g/2; Y:=g+Y; g:=F(T1,Y1)-F(T2,Y2); g:=(ih/ih1)*g/2; Y:=Y1+ih*(F(T1,Y1)+g+Y); g:=(ih1+ih2+ih3)/ih; g:=((ih1+ih2)/ih)*g; g:=(ih1/ih)*g/2; g:=(ih2/ih)*((ih2+ih3)/ih)/3+g; g:=(ih1/ih)*((3*ih1+4*ih2+2*ih3)/ih)/3+g; g:=(3*ih1+2*ih2+ih3)/(4*ih)+g; g:=(1/5+g)/24; Y:=Y+ih*ih*ih*ih*ih*g*PSI end {IntervalAB4}; procedure Find_stepsize (const p : Integer; const Y : array of interval; const lambda, PSI_width, Y_width, eps : Extended; const max_it : Integer; var h : array of Extended; out it, status : Integer); {------------------------------------------------------------------------------} { } { The function Find_stepsize calculates the step size for the next integration } { step for an interval method of Adams_Bashforth type of order p used to solve } { the initial value problem } { y' = 0.5y, y(0) = 1, } { with the exact solution y = exp(0.5t). To find this step the Newton } { iteration method is used. } { Data: } { p - order of method equal to the number k of method steps (p must } { be equal to 1, 2, 3 or 4), } { Y - array containing interval enclosures of solution at previous } { k steps (this array must be indexed from zero to k-1), } { lambda - constant occuring in the estimation } { width(F(T,Y)) <= lambda*(width(T)+width(Y)), } { PSI_width - width of PSI(delta_t, delta_y), } { Y_width - required width of Y at t_k, } { eps - accuracy for finding h for the next inegration step, } { max_it - maximal number of iterations, } { h - array (indexed from zero) containing step sizes at previous } { k-1 steps and an initial approximation of h for the next step } { included in h[k]. } { Results: } { h[k] - step size for the next integration step, } { it - number of iterations performed to find h[k], } { status - a variable which within the procedure is assigned the value of: } { 1 - if p < 1 or p > 4, } { 2 - if PSI_width <= 0 or any of h[l] for l = 1, 2, ... , k-1 is } { less than zero or Y_width < 0 or Y_width > 0.1 or } { eps <= 0 or eps > 0.1 } { 3 - if the required accuracy eps is not achieved in max_it } { iterations, } { 4 - if the calulated h[k] was less or equal to zero, } { 0 - otherwise. } { Note: If status = 1, 2, 3 or 4, then h[k] = 0 and it = 0. } { } {------------------------------------------------------------------------------} var k : Integer; h0, h1, q, q_prime, Yk_widths : Extended; begin status:=0; if (p<1) or (p>4) then status:=1 else if PSI_width<=0 then status:=2 else if (Y_width<0) or (Y_width>0.1) then status:=2 else if (eps<=0) or (eps>0.1) then status:=2 else begin k:=1; repeat if h[k]<0 then status:=2 else k:=k+1 until (k>p) or (status=2) end; if status<>0 then begin h[p]:=0; it:=0 end else begin Yk_widths:=0; for k:=p-1 downto 0 do Yk_widths:=Yk_widths+(k+1)*int_width(Y[k]); it:=0; h0:=h[p]; status:=3; repeat it:=it+1; case p of 1 : h1:=h0-(h0*h0*PSI_width/2+(1+h0*lambda)*Yk_widths-Y_width) /(h0*PSI_width+Y_width/2); 2 : begin if h0/h[p-1]<=1 then begin q:=h0; q_prime:=1; end else begin q:=h0*h0/h[p-1]; q_prime:=2*h0/h[p-1] end; h1:=h0-(h0*h0*(2*h0+3*h[p-1])*PSI_width/12 +lambda*q*Yk_widths+int_width(Y[p-1])-Y_width) /(h0*(h0+h[p-1])*PSI_width/2+lambda*q_prime*Yk_widths) end; 3 : begin if (h0/h[p-1]<=1) and (h0*(h0+h[p-1])/(h[p-1]*h[p-2])<=1) then begin q:=h0; q_prime:=1 end else if h0*(h0+h[p-1])/(h[p-1]*h[p-2])<=h0/h[p-1] then begin q:=h0*h0/h[p-1]; q_prime:=2*h0/h[p-1] end else begin q:=h0*h0*(h0+h[p-1])/(h[p-1]*h[p-2]); q_prime:=h0*(3*h0+2*h[p-1])/(h[p-1]*h[p-2]) end; h1:=h0-(h0*h0*((3*h0+4*(2*h[p-1]+h[p-2]))*h0 +6*h[p-1]*(h[p-1]+h[p-2]))*PSI_width/72 +lambda*q*Yk_widths+int_width(Y[p-1])-Y_width) /(h0*(2*h0*h0+2*h0*(2*h[p-1]+h[p-2])+h[p-1] *(h[p-1]+h[p-2]))*PSI_width/12+lambda*q_prime*Yk_widths) end; 4 : begin if (h0/h[p-1]<=1) and (h0*(h0+h[p-1])/(h[p-1]*h[p-2])<=1) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2]*(h[p-2]+h[p-3]))<=1) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2])*h[p-2]*h[p-3])<=1) then begin q:=h0; q_prime:=1 end else if (h0*(h0+h[p-1])/(h[p-1]*h[p-2])<=h0/h[p-1]) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2]*(h[p-2]+h[p-3]))<=h0/h[p-1]) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2])*h[p-2]*h[p-3])<=h0/h[p-1]) then begin q:=h0*h0/h[p-1]; q_prime:=2*h0/h[p-1] end else if (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2]*(h[p-2]+h[p-3])) <=h0*(h0+h[p-1])/(h[p-1]*h[p-2])) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2])*h[p-2]*h[p-3]) <=h0*(h0+h[p-1])/(h[p-1]*h[p-2])) then begin q:=h0*h0*(h0+h[p-1])/(h[p-1]*h[p-2]); q_prime:=h0*(3*h0+2*h[p-1]) /(h[p-1]*h[p-2]) end else if h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2])*h[p-2]*h[p-3]) <=h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2]*(h[p-2]+h[p-3])) then begin q:=h0*h0*(h0+h[p-1]) *(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2] *(h[p-2]+h[p-3])); q_prime:=h0*(h0 *(4*h0+6*h[p-1] +3*h[p-2]) +2*h[p-1] *(h[p-1]+h[p-2])) /(h[p-1]*h[p-2] *(h[p-2]+h[p-3])) end else begin q:=h0*h0*(h0+h[p-1]) *(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2]) *h[p-2]*h[p-3]); q_prime:=h0*(h0 *(4*h0+6*h[p-1] +3*h[p-2]) +2*h[p-1] *(h[p-1]+h[p-2])) /((h[p-1]+h[p-2]) *h[p-2]*h[p-3]) end; h1:=h0-(h0*h0*(((12*h0+15*(3*h[p-1]+2*h[p-2]+h[p-3]))*h0 +20*(h[p-1]*(h[p-1]+h[p-2])+(2*h[p-1]+h[p-2]) *(h[p-1]+h[p-2]+h[p-3])))*h0 +30*h[p-1]*(h[p-1]+h[p-2])*(h[p-1]+h[p-2]+h[p-3])) *PSI_width/1440 +lambda*q*Yk_widths+int_width(Y[p-1])-Y_width) /(h0*(12*h0*h0*h0+9*h0*h0*(3*h[p-1]+2*h[p-2]+h[p-3]) +8*h0*(h[p-1]*(h[p-1]+h[p-2]) +(2*h[p-1]+h[p-2])*(h[p-1]+h[p-2]+h[p-3])) +6*h[p-1]*(h[p-1]+h[p-2])*(h[p-1]+h[p-2]+h[p-3])) *PSI_width/288+lambda*q_prime*Yk_widths) end; end; if h1<=0 then status:=4 else if Abs(h1-h0)>eps then h0:=h1 else status:=0 until (it=max_it) or (status=0) or (status=4); if status=0 then h[p]:=h1 else begin h[p]:=0; it:=0 end end end {Find_stepsize}; var it, k, max_it, p, status : Integer; a, b_left, b_right, eps, lambda, PSI_width, real_h, real_t, w, Y_width : Extended; file_name, left, right, time : string; agree, output : Char; finish : Boolean; h : array of Extended; delta_t, delta_y, ih, ih1, ih2, ih3, T1, T2, T3, T4, Yh, Y1, Y2, Y3, Y4 : interval; Yk : array of interval; results : TextFile; time1, time2 : TDateTime; begin Writeln ('Solving the initial value problem'); Writeln ('y'' = 0.5y, y(0) = 1'); Writeln ('by interval k-step Adams-Bashforth methods of order p = k (p = 1, ' +'2, 3, 4)'); Writeln ('with variable step sizes for t in [0, a] and y in [b_left,' +' b_right]'); Writeln; Writeln ('Enter data:'); repeat Write (' order p = '); Readln (p); if (p<1) or (p>4) then Writeln ('p must be equal to 1, 2, 3 or 4') until (p>0) and (p<5); repeat Write (' a = '); Readln (a); if a<=0 then Writeln ('a must be greater then 0') until a>0; Write (' b_left = '); Readln (b_left); repeat Write (' b_right = '); Readln (b_right); if b_right<=b_left then Writeln ('b_right must be greater then b_left') until b_right>b_left; Str (a:26, right); delta_t.a:=0; delta_t.b:=right_read(right); Str (b_left:26, left); delta_y.a:=left_read(left); Str (b_right:26, right); delta_y.b:=right_read(right); SetLength (h, p+1); SetLength (Yk, p); for k:=1 to p do repeat Write (' h[', k, '] = '); Readln (h[k]); if h[k]<=0 then Writeln ('h[', k, '] must be greater then zero') until h[k]>0; repeat Write (' required Y width = '); Readln (Y_width); if (Y_width<=0) or (Y_width>0.1) then Writeln ('required Y width must be in the interval (0, 0.1]') until (Y_width>0) and (Y_width<=0.1); repeat Write ('accuracy for h, eps = '); Readln (eps); if (eps<=0) or (eps>0.1) then Writeln ('accuracy eps for finding h must be in the interval ' +'(0, 0.1]') until (eps>0) and (eps<=0.1); repeat Write (' max_it = '); Readln (max_it); if max_it<=0 then Writeln ('max number of iterations must be greater than 0') until max_it>0; lambda:=0.5; Writeln; repeat Write ('Select the output (s - screen, f - file) output = '); Readln (output) until (output='s') or (output='f'); if output='f' then begin repeat Writeln ('Enter the name of file'); Write ('(.txt will be added automatically) file_name = '); Readln (file_name); file_name:=file_name+'.txt'; Writeln ('The program will create a resulting file: ', file_name); Write ('Do you accept this name (y - yes, n - no) answer = '); Readln (agree) until agree='y'; Assign (results, file_name); Rewrite (results); Writeln (results, 'Solving the initial value problem'); Writeln (results, 'y'' = 0.5y, y(0) = 1'); Writeln (results, 'by an interval version of ', p, '-step ' +'Adams-Bashforth method of order p = ', p); Writeln (results, 'with variable step sizes for t in [0, a] and y ' +'in [b_left, b_right]'); Writeln (results, ' '); Writeln (results, 'a = ', a:4:2); Writeln (results, 'b_left = ', b_left:4:2); Writeln (results, 'b_right = ', b_right:4:2); for k:=1 to p do Writeln (results, 'h[', k:2, '] = ', h[k]:4:2); Writeln (results, 'required Y width = ', Y_width); Writeln (results, 'accuracy for h, eps = ', eps); Writeln (results, 'max_it = ', max_it) end; Writeln; Writeln ('Press Enter ...'); Readln; if output='s' then Write ('Enclosures of the exact solution obtained by the interval ' +'Adams-Bashforth method of ') else Write (results, 'Enclosures of the exact solution obtained by the ' +'interval Adams-Bashforth method of '); case p of 1 : begin if output='s' then begin Writeln ('first order'); Writeln end else begin Writeln (results, 'first order '); Writeln (results, ' ') end; k:=1; real_t:=0; T1:=0; iends_to_strings (T1, left, right); if output='s' then begin Writeln ('Initial intervals:'); Writeln ('T[0] = [', left, ',', right, ']') end else begin Writeln (results, 'Initial intervals:'); Writeln (results, 'T[0] = [', left, ',', right, ']') end; Y1:=1; Yk[0]:=Y1; iends_to_strings (Y1, left, right); if output='s' then begin Writeln ('Y[0] = [', left, ',', right, ']'); Writeln end else begin Writeln (results, 'Y[0] = [', left, ',', right, ']'); Writeln (results, ' ') end end; 2 : begin if output='s' then begin Writeln ('second order'); Writeln end else begin Writeln (results, 'second order '); Writeln (results, ' ') end; k:=2; real_t:=h[1]; Str (h[1]:26, left); ih1:=int_read(left); T1:=ih1; w:=Exp(0.5*real_t); Str (w:26, left); Y1:=int_read(left); Yk[1]:=Y1; T2:=0; Y2:=1; Yk[0]:=Y2; iends_to_strings (T2, left, right); if output='s' then begin Writeln ('Initial intervals:'); Writeln ('T[0] = [', left, ',', right, ']') end else begin Writeln (results, 'Initial intervals:'); Writeln (results, 'T[0] = [', left, ',', right, ']') end; iends_to_strings (Y2, left, right); if output='s' then Writeln ('Y[0] = [', left, ',', right, ']') else Writeln (results, 'Y[0] = [', left, ',', right, ']'); iends_to_strings (T1, left, right); if output='s' then Writeln ('T[1] = [', left, ',', right, ']') else Writeln (results, 'T[1] = [', left, ',', right, ']'); iends_to_strings (Y1, left, right); if output='s' then begin Writeln ('Y[1] = [', left, ',', right, ']'); Writeln end else begin Writeln (results, 'Y[1] = [', left, ',', right, ']'); Writeln (results, ' ') end end; 3 : begin if output='s' then begin Writeln ('third order'); Writeln end else begin Writeln (results, 'third order '); Writeln (results, ' ') end; k:=3; real_t:=h[1]+h[2]; Str (h[1]:26, left); ih1:=int_read(left); Str (h[2]:26, left); ih2:=int_read(left); Str (real_t:26, left); T1:=int_read(left); w:=Exp(0.5*real_t); Str (w:26, left); Y1:=int_read(left); Yk[2]:=Y1; T2:=ih1; w:=Exp(0.5*h[1]); Str (w:26, left); Y2:=int_read(left); Yk[1]:=Y2; T3:=0; Y3:=1; Yk[0]:=Y3; iends_to_strings (T3, left, right); if output='s' then begin Writeln ('Initial intervals:'); Writeln ('T[0] = [', left, ',', right, ']') end else begin Writeln (results, 'Initial intervals:'); Writeln (results, 'T[0] = [', left, ',', right, ']') end; iends_to_strings (Y3, left, right); if output='s' then Writeln ('Y[0] = [', left, ',', right, ']') else Writeln (results, 'Y[0] = [', left, ',', right, ']'); iends_to_strings (T2, left, right); if output='s' then Writeln ('T[1] = [', left, ',', right, ']') else Writeln (results, 'T[1] = [', left, ',', right, ']'); iends_to_strings (Y2, left, right); if output='s' then Writeln ('Y[1] = [', left, ',', right, ']') else Writeln (results, 'Y[1] = [', left, ',', right, ']'); iends_to_strings (T1, left, right); if output='s' then Writeln ('T[2] = [', left, ',', right, ']') else Writeln (results, 'T[2] = [', left, ',', right, ']'); iends_to_strings (Y1, left, right); if output='s' then begin Writeln ('Y[2] = [', left, ',', right, ']'); Writeln end else begin Writeln (results, 'Y[2] = [', left, ',', right, ']'); Writeln (results, ' ') end end; 4 : begin if output='s' then begin Writeln ('fourth order'); Writeln end else begin Writeln (results, 'fourth order '); Writeln (results, ' ') end; k:=4; real_t:=h[1]+h[2]+h[3]; Str (real_t:26, left); T1:=int_read(left); w:=Exp(0.5*real_t); Str (w:26, left); Y1:=int_read(left); Yk[3]:=Y1; w:=h[1]+h[2]; Str (w:26, left); T2:=int_read(left); w:=Exp(0.5*w); Str (w:26, left); Y2:=int_read(left); Yk[2]:=Y2; w:=h[1]; Str (w:26, left); T3:=int_read(left); w:=Exp(0.5*w); Str (w:26,left); Y3:=int_read(left); Yk[1]:=Y3; T4:=0; Y4:=1; Yk[0]:=Y4; iends_to_strings (T4, left, right); if output='s' then begin Writeln ('Initial intervals:'); Writeln ('T[0] = [', left, ',', right, ']') end else begin Writeln (results, 'Initial intervals:'); Writeln (results, 'T[0] = [', left, ',', right, ']') end; iends_to_strings (Y4, left, right); if output='s' then Writeln ('Y[0] = [', left, ',', right, ']') else Writeln (results, 'Y[0] = [', left, ',', right, ']'); iends_to_strings (T3, left, right); if output='s' then Writeln ('T[1] = [', left, ',', right, ']') else Writeln (results, 'T[1] = [', left, ',', right, ']'); iends_to_strings (Y3, left, right); if output='s' then Writeln ('Y[1] = [', left, ',', right, ']') else Writeln (results, 'Y[1] = [', left, ',', right, ']'); iends_to_strings (T2, left, right); if output='s' then Writeln ('T[2] = [', left, ',', right, ']') else Writeln (results, 'T[2] = [', left, ',', right, ']'); iends_to_strings (Y2, left, right); if output='s' then Writeln ('Y[2] = [', left, ',', right, ']') else Writeln (results, 'Y[2] = [', left, ',', right, ']'); iends_to_strings (T1, left, right); if output='s' then Writeln ('T[3] = [', left, ',', right, ']') else Writeln (results, 'T[3] = [', left, ',', right, ']'); iends_to_strings (Y1, left, right); if output='s' then begin Writeln ('Y[3] = [', left, ',', right, ']'); Writeln end else begin Writeln (results, 'Y[3] = [', left, ',', right, ']'); Writeln (results, ' ') end end end; Writeln; if output='f' then Writeln ('Calculating. Please, wait ...'); time1:=Now; case p of 1 : PSI_width:=int_width(PSI1(delta_t,delta_y)); 2 : PSI_width:=int_width(PSI2(delta_t,delta_y)); 3 : PSI_width:=int_width(PSI3(delta_t,delta_y)); 4 : PSI_width:=int_width(PSI4(delta_t,delta_y)) end; finish:=False; repeat Find_stepsize (p, Yk, lambda, PSI_width, Y_width, eps, max_it, h, it, status); if status=0 then begin real_h:=h[p]; if real_t+real_h