What solver/approach should I use for a system of pde's that pdepe seems unable to solve?

I'm trying to solve Poisson's equation, coupled to a non viscus Navier-Stokes & conservation of mass in 1D. V -> potential, v -> velocity, and n -> ion number density ne(V) is known.

% %%%%%%%%%%%%%%%%%%%%%% As the pde solver sees things %%%%%%%%%%%%%

%%% PDEs - v = u1, n = u2, V = u3

%  | 1 |     |u1|       |   0    |     | - u1*du1/dx - Ze/m * du3/dx | 
%  | 1 | * d |u2|  =  d |   0    |  +  |   - u1*du2/dx - u2*du1/dx   |
%  | 0 |   dt|u3|     dx| du3/dx |     |  - e/epsilon0( ne(V) - Zu2) |
%
%%%Boundary Conditions
%   
%  | u1 - ne(0) |  +  |0| .* |    0   |  =  | 0 |   neutral @ left boundary  
%  |     0      |  +  |1| .* | du3/dx |  =  | 0 |  Electric field = 0 @right

I recently noticed that there is a condition that one of the equations must be parabolic for pdepe. Which explains why the error: "This DAE appears to be of index greater than 1" in daeic12 in ode15s, is persistent. Regardless of what ic's I use and how self consistent they are.

I have been unable to find another appropriate solve/function. I am hoping not to have to write something completely from scratch as this is not my area of expertise.

It seems that the third eq. can be solved for each time step, and then the first two advanced explicitly, or with a grid that moves with v.

I am not looking forward to writing and developing this and am open to advice/suggestions.

Thanks