Если Вы имеете в виду, почему не решить получившуюся после дискретизации исходного уравнения систему

как

, то никто не запрещает так делать и это действительно будет быстрее. Преимуществами решения итерационным методом являются:
1) используется меньше оперативной памяти, поскольку не хранятся нулевые элементы матрицы

,
2) не надо использовать разреженные матрицы,
3) более наглядное представление разностного шаблона и задание граничных условий, чем при заполнении матрицы

и вектора

,
4) привычка решать данное уравнение именно этим методом.
Для ускорение работы функции выше можно использовать метод Гаусса-Зейделя, т.е. использовать формулу для phi(J,I):
Код:
phi(J,I) = 0.25*(phi_old(J+1,I)+phi_old(J,I+1)+phi(J-1,I)+phi(J,I-1)+rho(J,I)*dx^2/eps0);
или метод верхней релаксации, т.е. рассчитывать phi(J,I) как:
Код:
phi_tmp = 0.25*(phi_old(J+1,I)+phi_old(J,I+1)+phi(J-1,I)+phi(J,I-1)+rho(J,I)*dx^2/eps0);
phi(J,I) = omega*phi_tmp + (1-omega)*phi_old(J,I);
где

.