clear all factor = 1; n = 90; tval = [0 .3 .6 .9]; tfine = 0:.1:.9; h = 0.1; tvec(1) = 0; y(1) = 1; %y_approx = [1 1.029 1.131441 1.28742]; for j=1:9 tvec(j+1) = j*h; y(j+1) = y(j) + factor*h*(-y(j) + tvec(j) + 1); end y_approx(1) = y(1); y_approx(2) = y(4); y_approx(3) = y(7); y_approx(4) = y(10); coef = polyfit(tval,y_approx,2); poly = polyval(coef,tfine); p2 = coef(1); p1 = coef(2); p0 = coef(3); z(1) = p0; w(1) = 0; Tvec(1) = 0; Pvec(1) = 1; Dvec(1) = 0; for j=1:n dt = .9/n; t = (j-1)*dt; Tvec(j+1) = j*dt; P = p2*t^2 + p1*t + p0; D = 2*p2*t + p1 - (-P + t + 1); Dvec(j+1) = D; Pvec(j+1) = P; z(j+1) = z(j) + factor*dt*(-z(j) + t + 1 + D); w(j+1) = w(j) + factor*dt*(-w(j) + D); end max(abs(Pvec-z)) figure(1) plot(tval,y_approx,'or','linewidth',3) hold on plot(tfine,poly,'linewidth',3) plot(Tvec,z,'--g','linewidth',3) hold off set(gca,'Fontsize',[20]); xlabel('Time (s)') ylabel('Approximate and Pseudo-Problem Solution') legend('Approximate Solutions y_n','Polynomial Fit z(t)','Polynomial Approximations z_n','location','northwest') figure(2) plot(Tvec,w,'linewidth',3) hold on plot(Tvec,0*w,'linewidth',3) hold off set(gca,'Fontsize',[20]); xlabel('Time (s)') ylabel('Error w(t)')