It's good that you corrected that error replacing V with x. That's an important distinction in the integration process. However, it doesn't really change the difficult part of your problem, and that is the solution to your equation
p-q-lambda = sqrt(r*lambda+exp(lambda-x))
That is where you need to concentrate your energies.
As for the variable upper limit of integration, namely, a lambda which varies with x, both 'quad2d' and 'integral2' provide for the limits of integration on the inner integral to be variable in accordance with function handles. In other words they are able to integrate over a non-rectangular x,y region. You should read up the description of these function handles carefully to abide by their rules. That fact that your upper limit is not in the form of an explicit expression does not matter. If you can write a function that has the ability to return lambda's that correspond to entered x's, that is all that is needed to do the integration with either of these integration functions.
If I were doing your problem I would be tempted to use matlab's 'fzero' function to produce 'lambda'. Obviously you cannot do it with an anonymous function. It has to be an m-file function or subfunction. The integration routine is going to call on it with x values in the form of either vectors or perhaps matrices. You have to write the function to deal with each x value as a separate problem probably using for-loops and produce a corresponding lambda and then return with a like vector or matrix of these lambda values.
Within these for-loops you are faced with the same problem as before - given a single value of x, find the lambda value that satisfies the above equation. The 'fzero' routine allows you to enter a two-valued vector as your "estimate" argument, which I recommend. As the lower value you should use the least possible value that lambda could have. It should be so low that it always satisfies the inequality
p-q-lambda - sqrt(r*lambda+exp(lambda-x)) > 0
Perhaps zero would always serve there but you should check that out. It should also satisfy
r*lambda+exp(lambda-x) >= 0
so as to avoid complex results which would be disturbing to 'fzero'.
Choosing the upper value of the interval is perhaps a little trickier. The upper value must reverse the inequality above so that
p-q-lambda - sqrt(r*lambda+exp(lambda-x)) < 0
To accomplish this you may have to do a little iteration of your own, say repeatedly doubling some starting value until the inequality holds. Or perhaps you might discover some more sophisticated method. Just make sure it doesn't take any longer than the iteration done by 'fzero' does.
Once you have these two inequalities both satisfied, your solution will be pinned somewhere between them and 'fzero' will certainly "zero" in on the solution, and rather rapidly if I am not mistaken.
You will have to face the fact that computation will require more time than would be usual for a double integration procedure because of the heavy demands made on 'fsolve' to find 'lambda'.
That is an outline of how I would do the problem. I hope you manage to make it work.