% samplefd.m
% Mousumi Roy, Oct 2014
clear all
% finite difference solution of dy/dx = 2x^2 = f(x);  You all know what the
% answer should be!  lets use the boundary condition y=0.567 at x=0;

% lets define a domain and a discretization
dx = 0.01;        
x  = [dx:dx:1]'; % domain goes from 0 to 1
n  = length(x);  % there are n points

leftbc = 0.567;
% discretized system of equations for the case n=5;
% apply boundary condition at x=0; 
% y(0) = 0.567
% y(1) = y(0) + dx*f(x(1));
% y(2) = y(1) + dx*f(x(2));
% etc.
% y(5) = y(4) + dx*sin(x(5));
% or in matrix form:
%
% [ 1  0 0 0 0 ] [y(0)] = [leftbc      ]
% [ -1 1 0 0 0 ] [y(1)] = [dx*f(x(1))]
% [ 0 -1 1 0 0 ] [y(2)] = [dx*f(x(2))]
% [ 0  0 -1 1 0] [y(3)] = [dx*f(x(3))]
% [ 0  0 0 -1 1] [y(4)] = [dx*f(x(4))]

% let A be the matrix on the left made of ones and zeros
% let F be the matrix on the right
% we get A*y = F or, y = inv(A)*F

A = diag(ones(n,1)) - diag(ones(n-1,1), -1);
F = dx*2*x.^2;
% apply boundary condition
F(1,1) = leftbc;

y = inv(A)*F;

figure(1); clf;
plot(x,y,'ko'); hold on; plot(x,y,'k-');

figure(2); clf
%check
plot(x(1:n-1),diff(y)/dx, 'k-'); hold on
plot(x(2:n),F(2:n)/dx,'r-');