diff --git a/num/hw5/h-0.05.csv b/num/hw5/h-0.05.csv
diff --git a/num/hw5/h-0.1.csv b/num/hw5/h-0.1.csv
diff --git a/num/hw5/hw5.jl b/num/hw5/hw5.jl
@@ -1,15 +1,29 @@
+using PyPlot;
+
 asoln(x, y, t) = exp(0.75 * t) * sin(2 * x * y) * cosh(1.5 * (x + y));
 
+const debug_flag = true;
+
+macro if_debug(code)
+  if debug_flag
+    return code;
+  end
+end
+
 for h in [1/10; 1/20; 1/40]
-  const k = h;
+  w = open("h-$h.csv", "w");
+
+  const μ = 1.0;
+  const k = μ * h*h;
+  println("k = $k, h = $h");
 
   const aax = -k / (2 * h*h);
   const bbx = (k / (h*h) + 1);
   const ccx = aax;
 
   const aay = -k / (h*h);
   const bby = (2 * k / (h*h) + 1);
-  const ccy = aax;
+  const ccy = aay;
 
   xs = linspace(0.0, 1.0, Int(round(1.0 / h)));
   ys = copy(xs);
@@ -20,8 +34,22 @@ for h in [1/10; 1/20; 1/40]
   for (m, x) in zip(1:M, xs), (l, y) in zip(1:L, ys)
     u[m, l, 1] = asoln(x, y, 0);
   end
-
-  @time for n in 1:K-1
+  u_debug = (debug_flag) ? zeros(M, L, K) : zeros(1, 1, 1);
+  @if_debug(for (m, x) in zip(1:M, xs), (l, y) in zip(1:L, ys)
+    u_debug[m, l, 1] = asoln(x, y, 0);
+  end);
+
+  for n in 1:K-1
+    # calculate boundary terms
+    for l in 1:L
+      u[1, l, n+1] = asoln(0.0, ys[l], ts[n+1]);
+      u[M, l, n+1] = asoln(1.0, ys[l], ts[n+1]);
+    end
+    for m in 2:M-1
+      u[m, 1, n+1] = asoln(xs[m], 0.0, ts[n+1]);
+      u[m, L, n+1] = asoln(xs[m], 1.0, ts[n+1]);
+    end
+
     for l in 2:L-1 # consider boundary terms
       # calculate pi and qi for Thomas' algorithm
       p = zeros(L);
@@ -32,15 +60,15 @@ for h in [1/10; 1/20; 1/40]
                            (u[m, l+1, n] - 2 * u[m, l, n] + u[m, l-1, n]));
         denom = aax * p[m] + bbx;
         p[m+1] =  -ccx / denom;
-        q[m+1] = (ddx - aax * q[m]) / denom;
+        q[m+1] = (dd - aax * q[m]) / denom;
       end
       u[M, l, n+1] = asoln(1.0, ys[l], ts[n+1]);
-      for m=M-1:-1:1
+      for m=M-1:-1:2
         u[m, l, n+1] = p[m+1] * u[m+1, l, n+1] + q[m+1];
       end
     end
 
-    for m in 1:M
+    for m in 2:M-1
       # calculate pi and qi for Thomas' algorithm
       p = zeros(M);
       q = zeros(M);
@@ -51,18 +79,34 @@ for h in [1/10; 1/20; 1/40]
                          (u[m+1, l, n+1] - 2 * u[m, l, n+1] + u[m-1, l, n+1]));
         denom = aay * p[m] + bby;
         p[m+1] =  -ccy / denom;
-        q[m+1] = (ddy - aay * q[m]) / denom;
+        q[m+1] = (dd - aay * q[m]) / denom;
       end
       u[m, L, n+1] = asoln(xs[m], 1.0, ts[n+1]);
-      for m=M-1:-1:1
+      for l=L-1:-1:2
         u[m, l, n+1] = p[m+1] * u[m, l+1, n+1] + q[m+1];
       end
     end
 
-    err = zeros(M, L);
-    for (m, x) in zip(1:M, xs), (l, y) in zip(1:L, ys)
-      err[m, l] = u[m, l, n+1] - asoln(x, y, ts[n+1]);
+    @if_debug(begin
+        for l in 1:L
+          u_debug[1, l, n+1] = asoln(0.0, ys[l], ts[n+1]);
+          u_debug[M, l, n+1] = asoln(1.0, ys[l], ts[n+1]);
+        end
+        for m in 2:M-1
+          u_debug[m, 1, n+1] = asoln(xs[m], 0.0, ts[n+1]);
+          u_debug[m, L, n+1] = asoln(xs[m], 1.0, ts[n+1]);
+        end
+        for l in 2:L-1, m in 2:M-1
+          u_debug[m, l, n+1] = k / (h*h) * (u_debug[m+1, l, n] - 2 * u_debug[m, l, n] 
+                                            + u_debug[m-1, l, n] + 
+                                            2 * (u_debug[m, l+1, n] - 2 * u_debug[m, l, n] 
+                                                 + u_debug[m, l-1, n])) + u_debug[m, l, n];
+        end
+      end);
+
+    for m in 1:M, l in 1:L
+      write(w, "$(ts[n+1]), $(xs[m]), $(ys[l]), $(u[m, l, n+1]), $(asoln(xs[m], ys[l], ts[n+1])), $(u_debug[m, l, n+1])\n");
     end
-    println("Error at t=$(ts[n+1]): ", norm(err, 2));
   end
+  close(w);
 end