To the moderator: Please shift this to the appropriate forum if this isn't the place for it. EDIT -- It seems this (https://www.physicsforums.com/forumdisplay.php?f=189) is the right place for it? Sorry I overlooked this forum. ----- Hi everyone, I am solving some 2D boundary value problems in electromagnetics, using the Finite Element Method to familiarize myself with the method before I can apply it to my senior project. I have written a program in MATLAB which to solve the Laplace Equation using the FEM, on a rectangular solution domain defined by the coordinates (0,0), (1,0), (1,1), (0,1) [a square in this case]. The boundary conditions are V(x = 0) = V(x = 1) = 0 V(y = 0) = 0 V(y = 1) = V_0 (some fixed constant value) I have been able to write a program that generates a mesh for this region (please refer to attached diagram), and computes the value of the potential at every non-boundary node of the region (points of the triangular elements). Now, I need to use the element shape functions to interpolate the potential at points interior to an element using (a) the values of the potential at its 3 nodes (b) Lagrange polynomials (element shape functions) I want to obtain a contour plot for the solution in the entire solution domain. What is the best and easiest way to go about this? I thought of specifying several points within each triangular region and use the interpolation equation to calculate the values of the potential at those points. Then I would have a set of points (x, y, V(x,y)) which I could use to make a surface plot. I got a very poor plot when I tried this, so I am guessing there is something wrong in the way I've plotted this.