lambda = 10^9;
b = 2*pi*10^8;
D = 2*10^4;
F = 1.0;
R = 0.02;
beta = 10^-14;

alpha = D*beta/R;
gamma = -pi*F/b/R;

nx = 80;
ny = nx;

dx = lambda/nx;
dy = b/(ny);

xdist = [0:dx:lambda];
ydist = [0:dy:b];

[xg,yg] = meshgrid(xdist,ydist);

rhs = gamma*sin(pi*reshape(yg,(nx+1)*(ny+1),1)/b);
size(rhs)
rhs_gr = reshape(rhs,nx+1,ny+1);

laplace_operator = spdiags([((-2/dx^2 -2/dy^2)*ones((nx+1)*(ny+1),1)) ... 
                ((1/dy^2)*ones((nx+1)*(ny+1),1)) ...
                ((1/dy^2)*ones((nx+1)*(ny+1),1)) ...
                ((1/dx^2)*ones((nx+1)*(ny+1),1)) ...
                ((1/dx^2)*ones((nx+1)*(ny+1),1))], ...
		[0 -1 1 -(ny+1) (ny+1)],(nx+1)*(ny+1),(nx+1)*(ny+1));

% the x-derivative has only two diagonals (centered differences)
x_deriv_operator = spdiags([((1/2/dx)*ones((nx+1)*(ny+1),1)) ...
                     ((-1/2/dx)*ones((nx+1)*(ny+1),1))], ...
		     [(ny+1) -(ny+1)],(nx+1)*(ny+1),(nx+1)*(ny+1));

operator = laplace_operator + alpha*x_deriv_operator;

for i=2:ny
  rhs(i)   = 0;
  operator(i,i)   = 1;
  operator(i,i-1) = 0;
  operator(i,i+1) = 0;  
  operator(i,i+(ny+1)) = 0;    

  rhs(i+nx*(ny+1))   = 0;
  operator(i+nx*(ny+1),i+nx*(ny+1))   = 1;
  operator(i+nx*(ny+1),i+nx*(ny+1)-1) = 0;
  operator(i+nx*(ny+1),i+nx*(ny+1)+1) = 0;  
  operator(i+nx*(ny+1),i+nx*(ny+1)-(ny+1)) = 0;    
end

for i=2:nx
  rhs(1+(i-1)*(ny+1))   = 0;
  operator(1+(i-1)*(ny+1),1+(i-1)*(ny+1))   = 1;
  operator(1+(i-1)*(ny+1),1+(i-1)*(ny+1)-(ny+1)) = 0;
  operator(1+(i-1)*(ny+1),1+(i-1)*(ny+1)+(ny+1)) = 0;  
  operator(1+(i-1)*(ny+1),2+(i-1)*(ny+1)) = 0;    

  rhs((ny+1)+(i-1)*(ny+1))   = 0;
  operator((ny+1)+(i-1)*(ny+1),(ny+1)+(i-1)*(ny+1))   = 1;
  operator((ny+1)+(i-1)*(ny+1),(ny+1)+(i-1)*(ny+1)-(ny+1)) = 0;
  operator((ny+1)+(i-1)*(ny+1),(ny+1)+(i-1)*(ny+1)+(ny+1)) = 0;  
  operator((ny+1)+(i-1)*(ny+1),(ny)+(i-1)*(ny+1)) = 0;    
end

rhs(1)   = 0;
operator(1,1)      = 1;
operator(1,2)      = 0;
operator(1,ny+2)   = 0;

rhs(ny+1)   = 0;
operator(ny+1,ny+1)     = 1;
operator(ny+1,ny)       = 0;
operator(ny+1,2*(ny+1)) = 0;

rhs(1+(nx)*(ny+1))                = 0;
operator(1+(nx)*(ny+1),1+(nx)*(ny+1))   = 1;
operator(1+(nx)*(ny+1),2+(nx)*(ny+1))   = 0;
operator(1+(nx)*(ny+1),1+(nx-1)*(ny+1)) = 0;

rhs((nx+1)*(ny+1))                = 0;
operator((nx+1)*(ny+1),(nx+1)*(ny+1))   = 1;
operator((nx+1)*(ny+1),(nx)*(ny+1))     = 0;
operator((nx+1)*(ny+1),(nx+1)*(ny+1)-1) = 0;


psi = operator\rhs;

psi_gr = reshape(psi,nx+1,ny+1);
figure(2);contourf(xg,yg,psi_gr);title('solution?')
[v,u]=gradient(psi_gr,dx,dy);
u = -u;
figure(3)
ns = 2;
quiver(xg(1:ns:ny+1,1:ns:nx+1),yg(1:ns:ny+1,1:ns:nx+1), ...
	u(1:ns:ny+1,1:ns:nx+1),v(1:ns:ny+1,1:ns:nx+1),10)