Looks to be a very problem-specific solution with many (essentially all) BCs hardcoded as "magic numbers"; probably going to be somewhat difficult to generalize.
Only forces mentioned are initialized to zero; not clear, but I'm guessing that is solution, not input.
In the section
% converting uniform pressure into equivalent nodal forces. For this
% integration we are using one point Guass quadrature.
fe = Force(nnel,shape,P) ; % Force vector
there's a computation of a force vector from (one presumes a) function Force but there is no definition for Force in the posted code.
This appears as is to be written only for a uniform pressure load, P; you could approximate the solution by instead of using a uniform pressure over the full plate, turn P into a 2D array and set it to zero excepting at the element containing the desired point and setting P=F_applied/A_element for that element, zero elsewhere.
Where would be based on the fixed nodes as are defined at the beginning as to there only being 20x20 nodes. Of course, it wouldn't actually be a point load, but the approximation of a uniform pressure of same total force spread over the area bounding the desired location.
To actually include a point load the formulation would have to be modified significantly.
