function [e1,e2]=neumannbceg(N)

% solve heat equation on [0,1] with t>0, Dirichlet bc at x=1, u_x+u=0 at
% x=0, exact solution given by uexact function
% first use ficticious node in solution of mixed bc then use 2nd order
% method involving U0, U1 and U2

Tfinal=0.25;
theta=0.5;

% define and plot exact solution
uexact=@(x,t) exp(-(3*pi/2)^2*t).*(sin(3*pi*x/2)-3*pi/2*cos(3*pi*x/2));
figure(1), clf
fplot(@(x) uexact(x,Tfinal),[0,1],'b','LineWidth',2)
hold on

% set up grid in space and time
xx=linspace(0,1,N+1);
dx=xx(2);
mu=0.5;
M=N^2*Tfinal/mu;
tt=linspace(0,Tfinal,M+1);
%dt=tt(2);

% initial condition
U=uexact(xx',0);
g=zeros(N+1,1);

[LHSmat,RHSmat]=makemats(N,mu,theta,dx);

for m=2:M+1
    g(N+1)=uexact(1,tt(m));
    Unew=LHSmat\(RHSmat*U+g);
    U=Unew;
end
figure(1)
plot(xx,U,'r','LineWidth',2)
e1=max(abs(U-uexact(xx',0.25)));
figure(2), clf
plot(xx',abs(U-uexact(xx',0.25)),'r','LineWidth',2)
hold on

% now do second scheme for application of mixed bc
% initial condition
U=uexact(xx',0);
RHSmat(1,:)=0;
LHSmat(1,:)=0;
LHSmat(1,1)=2*dx-3;
LHSmat(1,2)=4;
LHSmat(1,3)=-1;

for m=2:M+1
    g(N+1)=uexact(1,tt(m));
    Unew=LHSmat\(RHSmat*U+g);
    U=Unew;
end
figure(1)
plot(xx,U,'g','LineWidth',2)
xlabel('x'),ylabel('u(x,0.25)')
legend('exact solution','numerical solution with ghost point','numerical solution','location','SouthEast')
e2=max(abs(U-uexact(xx',0.25)));
figure(2)
plot(xx',abs(U-uexact(xx',0.25)),'g','LineWidth',2)
xlabel('x'),ylabel('error in u(x,0.25)')
legend('error in numerical solution with ghost point','error in numerical solution')



function [LHSmat,RHSmat]=makemats(N,mu,theta,dx)

e=ones(N+1,1);
A=spdiags([e -2*e e],-1:1,N+1,N+1);
A(1,1)=-2*(1-dx); A(1,2)=2;
A(N+1,N)=0; A(N+1,N+1)=0;

LHSmat=speye(N+1)-mu*theta*A;
Ip=speye(N+1); Ip(N+1,N+1)=0;
RHSmat=Ip+mu*(1-theta)*A;