function fig224

%  plots exact solution elliptic equation

clear *
clf

ep=0.05; b=0;  a=1;

% set parameters
nx=40;
x=linspace(-1,1,nx);
y=x;

alpha=1;  beta=1;
chi=0.5*sqrt(alpha^2 + beta^2 - 4*ep)/ep;

nn=40;
pi2=2*pi;
b0 = (b-a)*ifc(0,pi2,0,ep)/pi2;
bess0 = besseli(0,chi);
for in=1:nn
    ap(in)=(b-a)*ifs(0,pi2,in,ep)/pi;
    bp(in)=(b-a)*ifc(0,pi2,in,ep)/pi;
    bess(in)=besseli(in,chi);
end;

U=zeros(nx,nx);
for ix=1:nx
    for iy=1:nx
        r=sqrt(x(ix)^2+y(iy)^2);
        if r < 1
            phi=atan2(y(iy),x(ix));
            rp=chi*r;
            s=b0*besseli(0,rp)/bess0;
            for in=1:nn
                t1=besseli(in,rp)*(ap(in)*sin(in*phi)+bp(in)*cos(in*phi))/bess(in);
                s=s+t1;
            end;
            u(ix,iy)=a+s*exp(-0.5*(x(ix)+y(iy))/ep);
        end;
    end;
end;

% get(gcf)
set(gcf,'position', [1788 1052 573 333]);

uu=U+3;
hold on
grid on
box on
%axes('position',[.13  .43  .8  .53])

surf(x,y,uu')
%surfc(x,y,u')
view(38, 22);
colormap gray

% get(gca)
set(gca,'ztick',[1 2 3]);
set(gca,'zticklabel',{'-2';'-1';'0'})

set(gca,'xtick',[-1 0 1]);
set(gca,'ytick',[-1 0 1]);

%axis([-1 1 -1 1 -3 0])

contour(x,y,uu',9,'k')
xlabel('x-axis','fontsize',14,'fontweight','bold')
ylabel('y-axis','fontsize',14,'fontweight','bold')
%zlabel('solution','fontsize',14,'fontweight','bold')
% set the fontsize to 12 for the plot  '
set(gca,'fontsize',12); 
hold off

function f=fs(phi,n,ep)
f = exp(0.5*(cos(phi)+sin(phi))/ep)*sin(n*phi);

function f=fc(phi,n,ep)
f = exp(0.5*(cos(phi)+sin(phi))/ep)*cos(n*phi);

function s=ifs(a,b,n,ep);
%  simpson's method
N=max(100,100*n);
xd=linspace(a,b,N+1);
h=xd(2)-xd(1);
s=fs(xd(1),n,ep);
for j=2:2:N
	f2=fs(xd(j+1),n,ep);
	s=s+4*fs(xd(j),n,ep)+2*f2;
end;
s=h*(s-f2)/3;

function s=ifc(a,b,n,ep);
%  simpson's method
N=max(100,100*n);
xd=linspace(a,b,N+1);
h=xd(2)-xd(1);
s=fc(xd(1),n,ep);
for j=2:2:N
	f2=fc(xd(j+1),n,ep);
	s=s+4*fc(xd(j),n,ep)+2*f2;
end;
s=h*(s-f2)/3;