function fig9_21

% uses 1st order symmetric differences

%  int(0.5*T*u_x^2-wu,x=0..1)
%  w=-x^r

% n=number of internal grid points
%n=3
n=20
r=4;

x=linspace(0,1,n+2);
h=x(2)-x(1);

T=1;
alpha=h^2/T;
x=linspace(0,1,n+2);
h=x(2)-x(1);

A=zeros(n,n);
A(1,1)=2; A(1,2)=-1;
for ix=2:n-1
    A(ix,ix-1)=-1; A(ix,ix)=2; A(ix,ix+1)=-1;
end
A(n,n-1)=-1; A(n,n)=2;

cond(A,inf)

for i=1:n
    w=-x(i+1)^r;
    b(i)=alpha*w;
end

u=A\b';
u=[0 u' 0];

ne=100;
xe=linspace(0,1,ne);
for i=1:ne
    ue(i)=xe(i)*(xe(i)^(r+1)-1)/((r+1)*(r+2)*T);
end

clf
% get(gcf)
set(gcf,'Position', [22 1107 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.06,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

plot(x,u,'o','MarkerSize',8,'LineWidth',2)
%plot(x,u,'--','MarkerSize',8,'LineWidth',2)
hold on
plot(xe,ue,'r','lineWidth',1.6)
grid on
set(gca,'ytick',[-0.02 -0.01 0]);
%set(gca,'YTickLabel',{'0.02';'0.01';'0.0'})
xlabel('x-axis')
ylabel('u-axis')
set(gca,'FontSize',16,'FontWeight','bold')

hold off

%exportgraphics(gcf,'/Users/mark/Desktop/beam2.eps')






















