%
% FDS.M Example of Finite Difference Method for Solving ODEs
%
% This demo uses the finite difference technique as the solution method
% for a particular ODE. The mesh interval is constant.
%
% Description of given differential equation
% y'' + 3xy' + 7y = cos(2x) with y(0) = 1 and y(pi) = 0
%
% A related Matlab file, FD.M, is also available (this contains a more detailed
% version of this example - compares two different approximations for y' and
% performs a mesh sensitivity study).
%
% File prepared by J. R. White, UMass-Lowell (Aug. 2003).
%
%
% getting started
clear all, close all, nfig = 0;
%
% set domain limits and boundary conditions
xo = 0; xf = pi; yxo = 1; yxf = 0;
%
% user defined number of unknowns in problem
N = input('Input number of unknowns in problem (N): ');
%
% compute interval size and discrete x vector
dx = (xf-xo)/(N+1); dx2 = dx*dx; x = (xo+dx):dx:(xf-dx);
%
%
% Approximations: y'' = [y(i-1) - 2y(i) + y(i+1)]/(dx^2) error => (dx^2)
% (Case A) y' = [y(i+1) - y(i-1)]/(2dx) error => (dx^2)
%
%
% setup matrix eqns. (treat boundary terms as special cases)
a = zeros(N,N); b = zeros(N,1);
for i = 2:N-1
a(i,i-1) = 1-3*x(i)*dx/2;
a(i,i) = -2+7*dx2;
a(i,i+1) = 1+3*x(i)*dx/2;
b(i) = dx2*cos(2*x(i));
end
% left boundary
a(1,1) = -2+7*dx2;
a(1,2) = 1+3*x(1)*dx/2;
b(1) = dx2*cos(2*x(1))-(1-3*x(1)*dx/2)*yxo;
% right boundary
a(N,N-1) = 1-3*x(N)*dx/2;
a(N,N) = -2+7*dx2;
b(N) = dx2*cos(2*x(N))-(1+3*x(N)*dx/2)*yxf;
% solve system of eqns
y = a\b;
% add boundary points to solution for plotting
za = [yxo y' yxf]; xa = [xo x xf];
%
% plot results
nfig = nfig+1; figure(nfig)
plot(xa,za,'r-','LineWidth',2)
axis([0 3.2 -2 1]);
title(['FDs: FD Method for Case A (',num2str(N),' pts)'])
xlabel('x values'),ylabel('y values'),grid
%
% end of demo