FEM.m code is below that needs to be altered as per part b:

% This code solves the BVP u''(x) + p(x) = 0 using Galerkin FEM clear all;

%%%%%%%% Preliminary variables %%%%%%%% N=5; % Number of node points (not elements). L=1; dx=L/(N-1); x=0:dx:L; p=ones(N,1).*x'; % This is p(x) in the ODE. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Create the trial functions psi and their derivatives dpsi psi=zeros(N,length(x)); dpsi=zeros(N,length(x)); [psi,dpsi] = psiterms(x, dx, N, psi);

% Having a look at the psi test functions figure(1); for i=1:N plot(x,psi(i,:),'color',rand(1,3)) hold on end axis([0 1 0 1]) xlabel('$x$','Interpreter','latex','Fontsize',28) ylabel('$\psi_{i}$','Interpreter','latex','Fontsize',28) title('Linear test functions $\psi_{i}$','Interpreter','latex','Fontsize',28) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculating the entries of the matrix A by simple integration % The diagonal entries A=zeros(N,N); for i=1:N if i==1 A(i,i)=dx*max(abs(dpsi(i,:)))^2.0; elseif i==N A(i,i)=dx*max(abs(dpsi(i,:)))^2.0; else A(i,i)=2.0*dx*max(abs(dpsi(i,:)))^2.0; end end % The off-diagonal entries for i=1:N for j=1:N if i~=j & abs(i-j)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculating the RHS vector f by simple trapezoidal integration of p*psi for i=1:N term(i,:)=p(:)'.*psi(i,:); % term is p(x)*psi(x). end

% Integrating f=zeros(N,1); for i=1:N for j = 2:N f(i) = f(i) + +(term(i,j)+term(i,j-1))/2*(x(j)-x(j-1)); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% We must delete the first and last rows and columns of A, the first % and last entries of f and the first and last rows of psi: %....

% Solving the linear system %....

% Calculate the full solution: %....

%Plot a comparison between exact and numerical solutions:

Problem 1 (The only problem!) The Matlab script FEM-m' has been created for you to have a play with. Make sure you also have the function file psiterms.m' in the same folder that you are running FEM.m from! Most of the hard work has been done for you already! It solves the BVP described in the notes, d2u (0) u(L) (1) dac2 for the simple choice of p(z) r, for 0 z S 1. The code has been set up with an initial choice of 5 rid elements that is 5 nodes at the points zi 0, 0.25, 0.5, 0.75 and 1 Therefore in this case N (A is a 5 x 5 matrix, f and u are vectors of length 5 x 1) and the code solves the linear system, where u is the vector of coefficients wi which are later used to build the FEM numerical solution in, as described in the background material. The code works by firstly defin- ing the r grid, step size and number of elements. It then calculates the linear trial functions whi, and their derivatives, as defined in the notes. The entries of the matrix A and the vector f are then determined by equations (4) and (5). Your task is as follows: Problem 1 (The only problem!) The Matlab script FEM-m' has been created for you to have a play with. Make sure you also have the function file psiterms.m' in the same folder that you are running FEM.m from! Most of the hard work has been done for you already! It solves the BVP described in the notes, d2u (0) u(L) (1) dac2 for the simple choice of p(z) r, for 0 z S 1. The code has been set up with an initial choice of 5 rid elements that is 5 nodes at the points zi 0, 0.25, 0.5, 0.75 and 1 Therefore in this case N (A is a 5 x 5 matrix, f and u are vectors of length 5 x 1) and the code solves the linear system, where u is the vector of coefficients wi which are later used to build the FEM numerical solution in, as described in the background material. The code works by firstly defin- ing the r grid, step size and number of elements. It then calculates the linear trial functions whi, and their derivatives, as defined in the notes. The entries of the matrix A and the vector f are then determined by equations (4) and (5). Your task is as follows