A new construction of fractal interpolation surfaces, using solutions of partial differential equations, is presented. We consider a set of interpolation points placed on a rectangular grid and a specific PDE, such that its Dirichlet's problem is uniquely solvable inside any given orthogonal region. We solve the PDE, using numerical methods, for a number of regions, to construct two functions H and B, which are then used to produce the fractal surface, as the attractor of an appropriately chosen recurrent iterated function system.