If you're in an academic setting, you can probably use the Silvaco/Simucad BSIM3/4 Verilog-A code ("not for commercial use" and you have to register).  Since you don't actually need the derivatives for a Newton solver (just evaluating Ids for given W,L -- and biases), you can probably pretty quickly convert the V-A code to MATLAB (MUCH more quickly than convering C code).   Prof. Roychowdhury at U Minn was using a Verilog-A version of BSIM3 to interface with some MATLAB functions for obtaining dc convergence using some new homotopy algorithms; you might try contacting him to see if he converted V-A to MATLAB or has some other way of evaluating V-A from MATLAB.

Hi, Peter -
BSIM3 has been used down even to 90nm CMOS, I believe.  Of course, depending on your application, you may find problems even in larger geometries; my company (Analog Devices) needed to enhance the model for high-precision circuits in 180nm (0.18u).

The code has a parameter TYPE that you can set to -1 to simulate PMOS.  And yes, you have to set (extract) a bunch of parameters for the model to match data for a given technology.

A fair number of the parameters are used for the diodes (body-drain, body-source), current and capacitance.  Oh, yes, and there's a lot of the code used for the CAPMOD channel charge/capacitance.  I see a comment in the copy I have
 // C-V model
(capacitance as a function of voltage), and you can delete that section (all the way down to and including the four lines after // Add depletion capacitance) if you just want Ids.  There aren't that many parameters used only in the C-V model, though, so you won't save much in that regard.

-Geoffrey