Thank you both @acer @Rouben Rostamian for your replies. 

I will repost. 

@Saha 

In your second problem you have not specified n, S and Nc. If you do that from what you say in your "a" defintions it runs for me. 

Always make sure there are no unkown variables in your boundary conditions otherwise it will not run. 

See attached. 

SCM2_response.mw

If you can upload your worksheet with your attempted solution, or provide more details to what you exactly are looking to solve this would be helpful. 

The big green arrow is used for uploading worksheets or the contents of your sheet. 

@filipm 

You have some issues in your physics. The final solution in your image is not the solution of a spacetime with that stress energy distribution. That solution solves the vacuum field equations with a cosmological constant. 

Please keep on eye on this thread as I will write an updated solution shortly. 

@ecterrab 

Thanks for reminding me of that option, I fail to remember that option exists when I do calculations with Lambda. 

Can you give some more information about the stress-energy tensor? It is incomplete currently but I believe you are using a perfect fluid? 

If you could also upload your attempts that would be great. 

Lastly here is the method for solving just the 2+1 with a cosmological constant and no stress-energy. Once you give more information of confirm it is a perfect fluid I will update the response. 

Response.mw

@sofreevique 

The green arrow that appears at the top of the text box (see image below) 

Click on it, go to choose files as it says and click on your saved worksheet then click upload and insert link. 

You will get the best responses if you post your worksheet to show what you have attempted. 

Or at the vary least post the ODEs you are working with. 

You can use the green arrow to upload your worksheet or images. 

@ecterrab 

Hello, I am sorry for replying so late. 

The desired output should be:

Details can be found in: https://physics.stackexchange.com/questions/93157/variation-of-modified-einstein-hilbert-action where some terms have been ignored. 

Thank you for taking a look into this. 

Just out of curiousity what are the exact ODEs you are trying to solve? 

@acer 

It does do rather interesting things in the negative region..something I had not thought about.

Is there any portion of the code that does deal with stepsize/speed? Not that I want to "speed up" the process just more of a curiousity on how it works. 

The positive result is also interesting, not what I expected. I will keep tinkering with the parameter ranges and see if anything can occur. I appreciate this worksheet, it appears my attempt was far too simplistic. 

Thanks!

@acer 

As mentioned in my comment above I am still have issues with CW suggestion and still recieving errors. This is a much more involed solution then I was expecting although much more elegant than my attempts. I have some questions however, some may be minor in nature. 

This method replaces my looping attempt by plotting the region of which the inequlaity you is satisfied in blue by examining the region alpha=-1.5..0 and beta -2..0 in step sizes of epsilon =e^-5? Is that the correct understanding?

This reason I ask if because I am interested in regions of positive alpha and beta only and switching to alpha=0..1.5 and beta=0..2 and running it has yet to produce a result after approximately 15 minutes real time. So if i wanted to "speed up" could I make the epsilon larger at the cost of precision? 

Thanks

@Christian Wolinski 

Using this substitution I still recieve an error - for some reason I cannot upload worksheet contents only the worksheet link itself. It appears to not be able to actually compute the integration as I recieve the new error in the modified sheet below. 

NestedError_Mod.mw

@Axel Vogt 

In your language I suppose the zero. In general H is less then G and that is why I want to find when H>G but writiing it as you do then I would want to find when f would be zero. 

My H is

Int(alpha^(3/2)*exp(-1/2*erf(1/2*sqrt(-2*alpha)*t)^2/erf(1/2*sqrt(-2*alpha)*b)^2 - 1/2*alpha*t^2)*I/(Pi*erf(1/2*sqrt(-2*alpha)*b)), t = -infinity .. infinity)and G is

and G is

-alpha/(4*Pi*erf(sqrt(-2*alpha)*b/2)^2)

alpha is in (10^(-9),1) and b is in (1,10^3). My procedure is found in LogLogProper.mw 

I am happy with what it returns and it does it quick, I just don't know if it is the best approach. 

If there is a better way to do it and you want to ammend it I appreciate it. Thank you. 

@acer 

I did not believe the end game was entirely prevelant for this question, I was originally just interested in the reason for the numeric integration issue but things spiraled from there. 

I also want to avoid the "do my work for me" reputation on this forum, all though I understand that this may be hard for a new user.