function reconstructed_img = pde_solve(pic)

rows = size(pic,1);
cols = size(pic,2);

disp('Generating Laplacian operator matrix');
T0 = clock;
A = Laplacian_mat(rows, cols);
etime(clock, T0)

disp('Modifying Laplacian operator matrix for boundary conditions');
T1 = clock;
b = reshape(pic', rows*cols, 1);
A = lap_bound_mat(A, b, rows, cols);
etime(clock, T1)

disp('Solving linear system');
T2 = clock;
reconstructed_img = reshape(A\b, cols, rows)';
%[L, U] = luinc(A, 1e-2);
%reconstructed_img = bicgstab(A, b, 1e-3, 100, L, U, b);
etime(clock, T2)