%  This program solves the fourth order two-point b.v.p.
%  
%         [au'']'' = f  in [0,1]
%  
%  with boundary conditions
%
%          u = u0,  u' = u1   (clamped)      
%  or      
%          u'' = 0,  (au'')' = 0  (free)
%
%  at  x = 0  and  x = 1,  respectively.
%  
%  using  C1  piecewice qubic elements and mesh size  h.
%  
%  
echo off
clear
clg
disp('This program solves the fourth order two-point b.v.p.')
disp('')
disp('       DD(aDDu) = f  in [0,1]')
disp('')
disp('with boundary condtions')
disp('')
disp('       u = u0,   Du = u1     (clamped)')
disp('or')
disp('       DDu = 0,    D(aDDu) = 0    (free)')
disp('')
disp('at  x=0  and  x=1,  respectively.')
disp('')
disp('using  C1  piecewice qubic elements and mesh size h.')
disp('')
disp('Give coefficient  a  in the equation')
acopy=input('a = ','s');
disp('Give rigth hand side  f  in the equation')
fcopy=input('f = ','s');
bctype=[1 1];
input('Boundary cond. type at  x=0,  c  (clamped)  or  f  (free) ?','s')
if ans=='c'; bctype(1)=0; u0(1)=input('u0 = '); u1(1)=input('u1 = '); 
elseif ans=='f'; bctype(1)=1;
end
input('Boundary cond. type at  x=1,  c  or  f ? ','s')
if ans=='c'; bctype(2)=0; u0(2)=input('u0 = '); u1(2)=input('u1 = '); 
elseif ans=='f'; bctype(2)=1;
end
% Ask for the exact solution, if given then plot it.
disp('Is the exact solution known?');
uexactflag=(input('answer y or n ','s') == 'y');
xxplot=0:.02:1;
[i,pltnr]=size(xxplot);
if uexactflag
  uexact =input('Give exact solution u(x) = ','s');
  Duexact=input('Give Du = du/dx = ','s');
  DDuexact=input('Give DDu = ','s');
  disp('Plotting the exact solution');
  x=xxplot(1);
  uplot(1)=eval(uexact);
  for i=2:pltnr
    x=xxplot(i); uplot(i)=eval(uexact);
  end
  hold off;  plot(xxplot,uplot,':');  hold on;
  title(['u(x) = ' uexact ' dotted; uh solid']);
end;
refine=1;       %Boolean used for refinement loop
first=1;        %Boolean used to indicate if turn in first loop
while refine
%Clear old data structures (otherwise, selecting next h > previous h 
% will cause an error)
clear h hh ni nn ndf stiff rhs f a u ucomp Ducomp DDucomp uex Duex DDuex err Derr DDerr xplot;
disp('Give desired mesh size h')
h=input('h = ');
ni=round(1/h);  h = 1/ni; %Modify h according to ni
disp(['Number of sub-intervals = ',int2str(ni)])
nn=(ni+1);
disp(['Number of nodes = ',int2str(nn)])
ndf=2*nn;
disp(['Number of degrees of freedom = ',int2str(ndf-2*sum(bctype))])
%
%disp('Before you use this progam, enter formulas for the four basis functions and their second derivatives on the unit interval [0,1]');
    b1copy='1-3*x^2+2*x^3';
    b2copy='x-2*x^2+x^3';
    b3copy='1-3*(1-x)^2+2*(1-x)^3';
    b4copy='-(1-x)+2*(1-x)^2-(1-x)^3';
    Db1copy='-6*x+6*x^2';
    Db2copy='1-4*x+3*x^2';
    Db3copy='6*(1-x)-6*(1-x)^2';
    Db4copy='1-4*(1-x)+3*(1-x)^2';
    DDb1copy='-6+12*x';
    DDb2copy='-4+6*x';
    DDb3copy='-6+12*(1-x)';
    DDb4copy='4-6*(1-x)';
%Set up six-diagonal stiffnes matrix
%
gauss=sqrt(3/5);
g(1)=1/2-gauss/2;
g(2)=1/2;
g(3)=1/2+gauss/2;
w(1)=5/18;
w(2)=8/18;
w(3)=5/18;
%
diag=zeros(size(1:ndf));
subdiag=zeros(size(1:ndf-1));
sub2diag=zeros(size(1:ndf-2));
sub3diag=zeros(size(1:ndf-3));
rhs=zeros(size(1:ndf));
%
for gp=1:3;
%
x=g(gp);
%
b1=eval(b1copy);
b2=eval(b2copy);
b3=eval(b3copy);
b4=eval(b4copy);
%
DDb1=eval(DDb1copy)/h^2;
DDb2=eval(DDb2copy)/h^2;
DDb3=eval(DDb3copy)/h^2;
DDb4=eval(DDb4copy)/h^2;
%
for i=1:ni;
x=(i-1)*h+g(gp)*h;
a(i)=eval(acopy);
f(i)=eval(fcopy);
end
%
hh=h*w(gp);
%
diag(3:2:ndf-1)=diag(3:2:ndf-1)+a.*DDb3.*DDb3.*hh;
diag(1:2:ndf-3)=diag(1:2:ndf-3)+a.*DDb1.*DDb1.*hh;
diag(4:2:ndf)=diag(4:2:ndf)+a.*DDb4.*DDb4.*hh;
diag(2:2:ndf-2)=diag(2:2:ndf-2)+a.*DDb2.*DDb2.*hh;
%
subdiag(3:2:ndf-1)=subdiag(3:2:ndf-1)+a.*DDb3.*DDb4.*hh;
subdiag(1:2:ndf-3)=subdiag(1:2:ndf-3)+a.*DDb1.*DDb2.*hh;
subdiag(2:2:ndf-2)=subdiag(2:2:ndf-2)+a.*DDb2.*DDb3.*hh;
%
sub2diag(1:2:ndf-3)=sub2diag(1:2:ndf-3)+a.*DDb1.*DDb3.*hh;
sub2diag(2:2:ndf-2)=sub2diag(2:2:ndf-2)+a.*DDb2.*DDb4.*hh;
%
sub3diag(1:2:ndf-3)=sub3diag(1:2:ndf-3)+a.*DDb1.*DDb4.*hh;
%
rhs(3:2:ndf-1)=rhs(3:2:ndf-1)+f.*b3.*hh;
rhs(1:2:ndf-3)=rhs(1:2:ndf-3)+f.*b1.*hh;
rhs(4:2:ndf)=rhs(4:2:ndf)+f.*b4.*hh;
rhs(2:2:ndf-2)=rhs(2:2:ndf-2)+f.*b2.*hh;
%
end %gp
%
supdiag=subdiag;
sup2diag=sub2diag;
sup3diag=sub3diag;
%
supdiag=[0 supdiag];
sup2diag=[0 0 sup2diag];
sup3diag=[0 0 0 sup3diag];
subdiag=[subdiag 0];
sub2diag=[sub2diag 0 0];
sub3diag=[sub3diag 0 0 0];
%
stiff=spdiags([sub3diag' sub2diag' subdiag' diag' supdiag' sup2diag' sup3diag'],-3:3,ndf,ndf);
stiff=full(stiff);
%
rhs=rhs';
if bctype(2)==0;
rhs=rhs-stiff(:,ndf-1)*u0(2)-stiff(:,ndf)*u1(2)*h;
rhs=rhs(1:ndf-2);
stiff=stiff(1:ndf-2,1:ndf-2);
ndf=ndf-2;
end
if bctype(1)==0;
rhs=rhs-stiff(:,1)*u0(1)-stiff(:,2)*u1(1)*h;
rhs=rhs(3:ndf);
stiff=stiff(3:ndf,3:ndf);
ndf=ndf-2;
end
%
rhs;
u=stiff\rhs;
u=u';
if bctype(1)==0;
u=[u0(1) u1(1)*h u];
end
if bctype(2)==0;
u=[u u0(2) u1(2)*h];
end
disp('Values of  u  at nodal points')
disp(u(1:2:2*ni+1)')
%
  % Graphics preparations
  %
  for gp=1:3;
  x=g(gp);
  b1=eval(b1copy);
  b2=eval(b2copy);
  b3=eval(b3copy);
  b4=eval(b4copy);
  Db1=eval(Db1copy)/h;
  Db2=eval(Db2copy)/h;
  Db3=eval(Db3copy)/h;
  Db4=eval(Db4copy)/h;
  DDb1=eval(DDb1copy)/h^2;
  DDb2=eval(DDb2copy)/h^2;
  DDb3=eval(DDb3copy)/h^2;
  DDb4=eval(DDb4copy)/h^2;
  ucomp(gp:3:3*ni-3+gp)=u(1:2:2*ni-1)*b1+u(2:2:2*ni)*b2+u(3:2:2*ni+1)*b3+u(4:2:2*ni+2)*b4;
  Ducomp(gp:3:3*ni-3+gp)=u(1:2:2*ni-1)*Db1+u(2:2:2*ni)*Db2+u(3:2:2*ni+1)*Db3+u(4:2:2*ni+2)*Db4;
  DDucomp(gp:3:3*ni-3+gp)=u(1:2:2*ni-1)*DDb1+u(2:2:2*ni)*DDb2+u(3:2:2*ni+1)*DDb3+u(4:2:2*ni+2)*DDb4;
  end
ucomp=[u(1) ucomp u(2*ni+1)];
  rn=1:ni;
  rn=rn*h;
  xplot(3:3:3*ni)=rn-ones(size(rn))*g(1)*h;
  xplot(2:3:3*ni-1)=rn-ones(size(rn))*g(2)*h;
  xplot(1:3:3*ni-2)=rn-ones(size(rn))*g(3)*h;
  xplot=[0 xplot 1];
  %Erase old plots?
figure(gcf)
  if first
    clgwin=1;     %If first time clear window and plot uexact, uh
    first=~first;
  else
    disp('Do you want to clear the graph window?');
    clgwin = (input('answer y or n ','s')=='y');
  end
  if clgwin
    if uexactflag
      hold off;  plot(xxplot,uplot,':',xplot,ucomp,'-');  hold on;
      title(['u(x) = ' uexact ' dotted;     uh solid']);
    else
      hold off;  plot(xplot,ucomp,'-');  hold on;
      title('uh');
    end
  else
    plot(xplot,ucomp,'-');
  end
  %
  % Error calculation
  %
  if uexactflag
   xplot=xplot(2:3*ni+1);
        for i=1:3*ni;
        x=xplot(i);
        uex(i)=eval(uexact);
        Duex(i)=eval(Duexact);
	DDuex(i)=eval(DDuexact);
	end
        DDucomp
        DDuex
        ucomp=ucomp(2:3*ni+1);
        err=max(abs(uex-ucomp));
	Derr=max(abs(Duex-Ducomp));
	DDerr=max(abs(DDuex-DDucomp));
	err
	Derr
	DDerr
  end; %of error calculation
  disp('Do you want to refine the mesh?');
  refine=(input('Answer y or n ','s')=='y');
end; %of refinement while loop
hold off
disp('If you want a hardcopy of the plot, enter the command print');
