function [a,b,c,d,e] = fill_y(JM,KM,B,dy,dt,jmax,kmax,epse)

a = zeros(jmax,kmax,3,3);
b = a;
c = a;
d = a;
e = a;

hy = 0.5*dt/dy;
he = dt/dy*epse;
for n = 1:3
  for m=1:3
    b(JM,KM,n,m) = -B(n,m)*hy;
    d(JM,KM,n,m) = B(n,m)*hy;
    c(JM,KM,n,m) = 0.0;
  end
  b(JM,KM,n,n) = b(JM,KM,n,n) - he;
  d(JM,KM,n,n) = d(JM,KM,n,n) - he;
  c(JM,KM,n,n) = 1.0 + c(JM,KM,n,n) + 2*he;
end