% Fixed-Free 1D Finite Difference Model using n nodes and f=1
% David Bloore, July 2011

% This code does not compensate for the error in u.  Can you fix it?

clc;
clear all;

% here the user is queried for the desired number of mesh points

n=input('N?');

% construction of h(n) and the domain for u(xx)

h=1/(n+1);
xx=linspace(0,1,1000);

% construction of H (starting from K)

K=sparse(toeplitz([2 -1 zeros(1,n-1)]));
H=K;
H(n+1,n+1)=1;

% don't forget the h^2 in f!

f=h^2*ones(n+1,1);

% solution to system, both linear algebra and analytic

u=H\f;
yy=-(1/2)*xx.^2+xx;

% we actually have to take u and put it in a larger vector so the line
% representing the solution starts at the origin

x=linspace(0,1,n+2);
y=zeros(1,n+2);
y(2:n+2)=u;
ya=-(1/2)*x.^2+x;

% computaion of error

err=y-ya;

% plot and parameters

plot(x,y,'--',xx,yy,'-',x,y,'o',x,err,'-.',x,err,'o');
legend('Solution to Hu=f','f(x)=-(1/2)x^2+(1/2)x','Mesh Point Values','Error','Errors at Mesh Points','Location','NorthWest');
title(['Fixed-Free System with f(x)=1, solved using Hu=f with ' sprintf('%d',n+1) ' mesh points']);
xlabel('0 \leq x \leq 1');
ylabel('Displacement from Equilibruim Position');
grid on;