M.file or external function to solve the following model

If you are interested in replication you may also want to take a look at this: