Higher-order numerical techniques are developed for the solution of (i) homogeneous heat equation u t = u xx and (ii) inhomogeneous heat equation u t = u xx + s(x, t) subject to initial condition u(x, 0) = f (x), 0 < x < 1, boundary condition u(0, t) = g(t)0 < t ≤ T and with non-local boundary condition(s) (i) b 0 u(x, t)dx = M (t) 0 < t ≤ T, 0 < b < 1 (ii) u(0, t) = (iii) u(1, t) = 1 0 φ(x, t)u(x, t)dx + g 1 (t), 0 < t ≤ T and 1 0 ψ(x, t)u(x, t)dx + g 2 (t), 0 < t ≤ T as appropriate. The integral conditions are approximated using Simpson’s 1 3 rule while the space derivatives are approximated by higher-order finite difference approxi- mations. Then method of lines, semidiscritization approach, is used to trans- form the model partial differential equations into systems of first-order linear ordinary differential equations whose solutions satisfy recurrence relations in- volving matrix exponential functions. The methods are higher-order accurate in space and time and do not require the use of complex arithmetic. Parallel algorithms are also developed and implemented on several problems from lit- erature and are found to be highly accurate. Solutions of these problems are compared with the exact solutions and the solutions obtained by alternative techniques where available.