% Finite Difference Code for BVP with Pseudo-Problem clear all N = 10; deg = input('Polynomial degree = '); h = 3/N; xs = h:h:3-h; x = 0:h:3; u = zeros(N+1,1); z = zeros(N+1,1); v0 = (2+(h^2)*(pi^2))*ones(N-1,1); v1 = -ones(N-2,1); A = diag(v0,0) + diag(v1,1) + diag(v1,-1); f = (h^2)*(pi^2)*sin((pi/3)*xs'); us = A\f; u(2:N,:) = us; Pval = polyfit(x',u,deg); Ps = polyval(Pval,x); Pss = polyval(Pval,xs); figure(1) plot(x,u,x,Ps,'--r','linewidth',2) set(gca,'Fontsize',[20]); xlabel('Time (s)') legend('Approximate Solution','Polynomial Fit','location','South') % % Solution to Pseudo- or nearby problem % for k = 1:deg-1 coef(k) = (deg-(k-1))*(deg-k); P2der(k) = coef(k)*Pval(k); end Ps2der = polyval(P2der,xs); F = -(h^2)*Ps2der' + (h^2)*(pi^2)*Pss'; zs = A\F; z(2:N,:) = zs; max(abs(z-u)) figure(2) plot(x,u,x,z,'--','linewidth',2) set(gca,'Fontsize',[20]); xlabel('Time (s)') legend('Approximate Solution','Pseudoproblem Solution','location','South') figure(3) plot(xs,f,xs,F,'--r','linewidth',2) set(gca,'Fontsize',[20]); xlabel('Time (s)') legend('pi^2 h^2 sin(pi x/3)','-h^2 D^2 P(x) - h^2 pi^2 P(x)','location','South')